######################
# Replication Nationscape Survey Figures, Models, & Appendix Tables.
# "Embracing the Status Hierarchy: How Immigration Attitudes, Prejudice, and
# Sexism Shaped Non-White Support for Trump”
# Geiger and Reny 2024
# Perspectives in Politics
# For questions please contact Jessica Geiger
# jessica.geiger@cgu.edu
######################

rm(list=ls())

pacman::p_load(haven, tidyverse, lubridate, stargazer, estimatr, simcf)
select <- dplyr::select

# simcf is not on CRAN and needs to be downloaded  
# though analyses could be closely replicated using another package like
# marginaleffects on CRAN

#setwd('~/replication')
setwd('/replication/')

#read in data
df <- read_csv('/data/ns_replication.csv')

# custom theme for plots
goldenScatterCAtheme <- theme(
  ## Removes main plot gray background
  panel.background = element_rect(fill = "white"),
  ## Golden rectangle plotting area (leave out for square)
  aspect.ratio = ((1 + sqrt(5))/2)^(-1),
  ## All axes changes
  axis.ticks.length = unit(0.5, "char"), # longer ticks
  ## Horizontal axis changes
  axis.line.x.top = element_line(size = 0.2), # thinner axis lines
  axis.line.x.bottom = element_line(size = 0.2), # thinner axis lines
  axis.ticks.x = element_line(size = 0.2), # thinner ticks
  axis.text.x = element_text(color = "black", size = 12),
  ## match type of axis labels and titles
  axis.title.x = element_text(size = 12,
                              margin = margin(t = 7.5, r = 0, b = 0, l = 0)),
  ## match type; pad space between title and labels
  ## Vertical axis changes
  axis.ticks.y = element_blank(), # no y axis ticks (gridlines suffice)
  axis.text.y = element_text(color = "black", size = 12,
                             margin = margin(t = 0, r = -4, b = 0, l = 0)),
  ## match type of axis labels and titles, pad
  axis.title.y = element_text(size = 12,
                              margin = margin(t = 0, r = 7.5, b = 0, l = 0)),
  ## match type of axis labels and titles, pad
  ## Legend
  legend.key = element_rect(fill = NA, color = NA),
  ## Remove unhelpful gray background
  ## Gridlines (in this case, horizontal from left axis only
  panel.grid.major.x = element_blank(),
  panel.grid.major.y = element_line(color = "gray45", size = 0.2),
  ## Faceting (small multiples)
  strip.background = element_blank(),
  ## Remove unhelpful trellis-like shading of titles
  strip.text.x = element_text(size=12), # Larger facet titles
  strip.text.y = element_blank(), # No titles for rows of facets
  strip.placement = "outside", # Place titles outside plot
  panel.spacing.x = unit(1.25, "lines"), # Horizontal space b/w plots
  panel.spacing.y = unit(1, "lines") # Vertical space b/w plots
)
class(goldenScatterCAtheme)

######################
# Figure 1 - Support for Status Hierarchies Across Racial Groups. 
######################

immpol <- df %>%
  group_by(race_cat) %>%
  summarise(y=weighted.mean(imm_policy,na.rm=T,w=weight)) %>%
  na.omit() %>%
  mutate(outcome='Immigration Policy')
ethno <- df %>%
  group_by(race_cat) %>%
  summarise(y=weighted.mean(ethnocentrism2,na.rm=T,w=weight)) %>%
  na.omit() %>%
  mutate(outcome='White Ethnocentrism')
rrdat <- df %>%
  group_by(race_cat) %>%
  summarise(y=weighted.mean(rr,na.rm=T,w=weight)) %>%
  na.omit() %>%
  mutate(outcome='Racial Resentment')
sexism <- df %>%
  group_by(race_cat) %>%
  summarise(y=weighted.mean(sexism,na.rm=T,w=weight)) %>%
  na.omit() %>%
  mutate(outcome='Sexism')

bind_rows(immpol, rrdat, ethno, sexism) %>%
  mutate(race_cat = factor(race_cat, levels=rev(c(
    'Black','Latino','Asian','White'
  ))), outcome=factor(outcome, levels=c(
    'Immigration Policy',
    'Racial Resentment',
    'White Ethnocentrism',
    'Sexism'
  ))) %>%
  ggplot(aes(x=race_cat, y=y)) +
  geom_segment(aes(x=race_cat, xend=race_cat, y=0, yend=y)) +
  geom_point(shape=21, size=3, fill='white') +
  facet_wrap(~outcome, ncol=4,strip.position = 'top') +
  goldenScatterCAtheme +
  labs(x='', y='Scale Value (0-1)') +
  scale_y_continuous(limits=c(0,0.6), breaks=c(0,0.2,0.4,0.6)) +
  scale_fill_brewer(palette = 'Set1', name='')  +
  theme(legend.position='none')
ggsave(width=10, height=4,filename = 'figures/figure1.eps')

######################
# Figure 2 - Immigration Attitudes and Trump Support.
######################

set.seed(12345)
form1 <- trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5
form2 <- trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5
form3 <- trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + korean + aapi_other
form4 <- trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
dat1 <- df[df$white==1,]
dat2 <- df[df$black==1,]
dat3 <-  df[df$aapi==1,]
dat4 <- df[df$latino==1,]

summary(m1 <- glm(form1, dat1 , family=binomial(link='logit')))
summary(m2 <- glm(form2, dat2, family=binomial(link='logit')))
summary(m3 <- glm(form3, dat3, family=binomial(link='logit')))
summary(m4 <- glm(form4, dat4, family=binomial(link='logit')))

pe1 <- coef(m1); vc1 <- vcov(m1); sb1 <- MASS::mvrnorm(10000, pe1, vc1)
pe2 <- coef(m2); vc2 <- vcov(m2); sb2 <- MASS::mvrnorm(10000, pe2, vc2)
pe3 <- coef(m3); vc3 <- vcov(m3); sb3 <- MASS::mvrnorm(10000, pe3, vc3)
pe4 <- coef(m4); vc4 <- vcov(m4); sb4 <- MASS::mvrnorm(10000, pe4, vc4)

xhyp <- seq(0,1,length.out=100)
mm1 <- cfMake(form1, dat1, nscen=length(xhyp))
mm2 <- cfMake(form2, dat2, nscen=length(xhyp))
mm3 <- cfMake(form3, dat3, nscen=length(xhyp))
mm4 <- cfMake(form4, dat4, nscen=length(xhyp))

mm1 <- cfChange(mm1, 'imm_policy', x=xhyp, scen=1:length(xhyp))
mm2 <- cfChange(mm2, 'imm_policy', x=xhyp, scen=1:length(xhyp))
mm3 <- cfChange(mm3, 'imm_policy', x=xhyp, scen=1:length(xhyp))
mm4 <- cfChange(mm4, 'imm_policy', x=xhyp, scen=1:length(xhyp))

res1 <- as_tibble(logitsimev(mm1, sb1)) %>% mutate(x=xhyp, group='White')
res2 <- as_tibble(logitsimev(mm2, sb2)) %>% mutate(x=xhyp, group='Black')
res3 <- as_tibble(logitsimev(mm3, sb3)) %>% mutate(x=xhyp, group='Asian')
res4 <- as_tibble(logitsimev(mm4, sb4)) %>% mutate(x=xhyp, group='Latino')

dat <- bind_rows(res1, res2, res3, res4) %>%
  mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))
diff_imm <- dat %>%
  filter(x %in% c(0,1)) %>%
  select(pe, group) %>%
  mutate(z=rep(c('min','max'),4)) %>%
  pivot_wider(id_cols = group,names_from=z,values_from=pe) %>%
  mutate(diff=max-min, outcome='Immigration') %>%
  select(group, diff, outcome)
labelling <- dat %>% filter(x==0 | x==1) %>% mutate(pe = round(pe, 2)) %>%
  mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))
dat %>%
  ggplot(aes(x=x, y=pe, ymin=lower, ymax=upper)) +
  geom_ribbon(alpha=0.3, fill='purple') +
  geom_line() + 
  facet_wrap(~group, ncol=4) +
  labs(x='Immigration Policy Attitudes', y='Pr(Vote Trump)') +
  scale_x_continuous(limits=c(-0.05,1.05),
                     breaks=c(0,0.25,0.5,0.75,1),
                     labels=c('0','0.25','0.5','0.75','1')) +
  scale_y_continuous(breaks=c(0,0.22,0.4,0.6,0.8,1),
                     limits=c(-0.05,1)) +
  goldenScatterCAtheme +
  geom_point(data=labelling, aes(x=x, y=pe,  group=group), size=7.5, shape=21, fill='white') +
  geom_text(data=labelling, aes(x=x, y=pe,  group=group, label=pe*100),size=2.5) +
  theme(axis.text.y = element_blank())
ggsave(width=8, height=4,filename = 'figures/figure2.eps')

######################
# Figure 3 - Racial Resentment and Trump Support.
######################

set.seed(12345)
form1 <- trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5
form2 <- trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5
form3 <- trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + korean + aapi_other
form4 <- trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other

summary(m1 <- glm(form1, dat1 , family=binomial(link='logit')))
summary(m2 <- glm(form2, dat2, family=binomial(link='logit')))
summary(m3 <- glm(form3, dat3, family=binomial(link='logit')))
summary(m4 <- glm(form4, dat4, family=binomial(link='logit')))

pe1 <- coef(m1); vc1 <- vcov(m1); sb1 <- MASS::mvrnorm(10000, pe1, vc1)
pe2 <- coef(m2); vc2 <- vcov(m2); sb2 <- MASS::mvrnorm(10000, pe2, vc2)
pe3 <- coef(m3); vc3 <- vcov(m3); sb3 <- MASS::mvrnorm(10000, pe3, vc3)
pe4 <- coef(m4); vc4 <- vcov(m4); sb4 <- MASS::mvrnorm(10000, pe4, vc4)

xhyp <- seq(0,1,length.out=100)
mm1 <- cfMake(form1, dat1, nscen=length(xhyp))
mm2 <- cfMake(form2, dat2, nscen=length(xhyp))
mm3 <- cfMake(form3, dat3, nscen=length(xhyp))
mm4 <- cfMake(form4, dat4, nscen=length(xhyp))

mm1 <- cfChange(mm1, 'rr', x=xhyp, scen=1:length(xhyp))
mm2 <- cfChange(mm2, 'rr', x=xhyp, scen=1:length(xhyp))
mm3 <- cfChange(mm3, 'rr', x=xhyp, scen=1:length(xhyp))
mm4 <- cfChange(mm4, 'rr', x=xhyp, scen=1:length(xhyp))

res1 <- as_tibble(logitsimev(mm1, sb1)) %>% mutate(x=xhyp, group='White')
res2 <- as_tibble(logitsimev(mm2, sb2)) %>% mutate(x=xhyp, group='Black')
res3 <- as_tibble(logitsimev(mm3, sb3)) %>% mutate(x=xhyp, group='Asian')
res4 <- as_tibble(logitsimev(mm4, sb4)) %>% mutate(x=xhyp, group='Latino')

dat <- bind_rows(res1, res2, res3, res4) %>%
  mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))

diff_rr <- dat %>%
  filter(x %in% c(0,1)) %>%
  select(pe, group) %>%
  mutate(z=rep(c('min','max'),4)) %>%
  pivot_wider(id_cols = group,names_from=z,values_from=pe) %>%
  mutate(diff=max-min, outcome='Racial Resentment') %>%
  select(group, diff, outcome)

labelling <- dat %>% filter(x==0 | x==1) %>% mutate(pe = round(pe, 2)) %>%
  mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))
dat %>%
  ggplot(aes(x=x, y=pe, ymin=lower, ymax=upper)) +
  geom_ribbon(alpha=0.3, fill='purple') +
  geom_line() + 
  facet_wrap(~group, ncol=4) +
  labs(x='Racial Resentment', y='Pr(Vote Trump)') +
  scale_x_continuous(limits=c(-0.05,1.05),
                     breaks=c(0,0.25,0.5,0.75,1),
                     labels=c('0','0.25','0.5','0.75','1')) +
  scale_y_continuous(breaks=c(0,0.2,0.4,0.6,0.8,1),
                     limits=c(-0.05,1)) +
  goldenScatterCAtheme +
  geom_point(data=labelling, aes(x=x, y=pe,  group=group), size=7.5, shape=21, fill='white') +
  geom_text(data=labelling, aes(x=x, y=pe,  group=group, label=pe*100),size=2.5) +
  theme(axis.text.y = element_blank())
ggsave(width=8, height=4,filename = 'figures/figure3.eps')

######################
# Appendix C1 - Trump Vote: Conservative Immigration Attitudes and Racial Resentment
# this code just saves the regressions from above in table format for the appendix
######################

summary(m1 <- glm(trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy +ideo5, data=df[df$white==1,],family=binomial(link='logit')))
summary(m2 <- glm(trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy +ideo5, data=df[df$black==1,],family=binomial(link='logit')))
summary(m3 <- glm(trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean, df[df$aapi==1,],family=binomial(link='logit')))
summary(m4 <- glm(trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy +ideo5+ cuban + puertorican + hispanic_other, data=df[df$latino==1,],family=binomial(link='logit')))
summary(m5 <- glm(trump_hr ~ rr +college_educ + household_income + female + age + pid7_legacy +ideo5, data=df[df$white==1,],family=binomial(link='logit')))
summary(m6 <- glm(trump_hr ~ rr +college_educ + household_income + female + age + pid7_legacy +ideo5, data=df[df$black==1,],family=binomial(link='logit')))
summary(m7 <- glm(trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean, df[df$aapi==1,],family=binomial(link='logit')))
summary(m8 <- glm(trump_hr ~ rr +college_educ + household_income + female + age + pid7_legacy +ideo5+cuban + puertorican + hispanic_other, data=df[df$latino==1,],family=binomial(link='logit')))

stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          out='tables/c1.tex',
          type="text",
          dep.var.labels = "",
          column.labels = c("White","Black", "AAPI","Latino", 
                            "White","Black", "AAPI","Latino",
                            "White","Black", "AAPI","Latino"),
          model.names = F,
          omit.stat = c('ll','aic'),
          title = "Logistic Regression Models Predicting Trump Vote",
          covariate.labels = c("Immigration Policy","Racial Resentment","College Education",
                               "Household Income","Gender", "Age", "Party ID", "Ideology",
                               "Indian", "Other AAPI", "Korean",
                               "Cuban", "Puerto Rican", "Hispanic Other", 
                               "Constant"))

######################
# Figure 4 - White Ethnocentrism and Trump Support.
######################

set.seed(12345)
form1 <- trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5
form2 <- trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5
form3 <- trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + korean + aapi_other
form4 <- trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other

summary(m1 <- glm(form1, dat1 , family=binomial(link='logit')))
summary(m2 <- glm(form2, dat2, family=binomial(link='logit')))
summary(m3 <- glm(form3, dat3, family=binomial(link='logit')))
summary(m4 <- glm(form4, dat4, family=binomial(link='logit')))

pe1 <- coef(m1); vc1 <- vcov(m1); sb1 <- MASS::mvrnorm(10000, pe1, vc1)
pe2 <- coef(m2); vc2 <- vcov(m2); sb2 <- MASS::mvrnorm(10000, pe2, vc2)
pe3 <- coef(m3); vc3 <- vcov(m3); sb3 <- MASS::mvrnorm(10000, pe3, vc3)
pe4 <- coef(m4); vc4 <- vcov(m4); sb4 <- MASS::mvrnorm(10000, pe4, vc4)

xhyp <- seq(0,1,length.out=100)
mm1 <- cfMake(form1, dat1, nscen=length(xhyp))
mm2 <- cfMake(form2, dat2, nscen=length(xhyp))
mm3 <- cfMake(form3, dat3, nscen=length(xhyp))
mm4 <- cfMake(form4, dat4, nscen=length(xhyp))

mm1 <- cfChange(mm1, 'ethnocentrism2', x=xhyp, scen=1:length(xhyp))
mm2 <- cfChange(mm2, 'ethnocentrism2', x=xhyp, scen=1:length(xhyp))
mm3 <- cfChange(mm3, 'ethnocentrism2', x=xhyp, scen=1:length(xhyp))
mm4 <- cfChange(mm4, 'ethnocentrism2', x=xhyp, scen=1:length(xhyp))

res1 <- as_tibble(logitsimev(mm1, sb1)) %>% mutate(x=xhyp, group='White')
res2 <- as_tibble(logitsimev(mm2, sb2)) %>% mutate(x=xhyp, group='Black')
res3 <- as_tibble(logitsimev(mm3, sb3)) %>% mutate(x=xhyp, group='Asian')
res4 <- as_tibble(logitsimev(mm4, sb4)) %>% mutate(x=xhyp, group='Latino')

dat <- bind_rows(res1, res2, res3, res4) %>%
  mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))

diff_ethno <- dat %>%
  filter(x %in% c(0,1)) %>%
  select(pe, group) %>%
  mutate(z=rep(c('min','max'),4)) %>%
  pivot_wider(id_cols = group,names_from=z,values_from=pe) %>%
  mutate(diff=max-min, outcome='White Ethnocentrism') %>%
  select(group, diff, outcome)

labelling <- dat %>% filter(x==0 | x==1) %>% mutate(pe = round(pe, 2)) %>%
  mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))
dat %>%
  ggplot(aes(x=x, y=pe, ymin=lower, ymax=upper)) +
  geom_ribbon(alpha=0.3, fill='purple') +
  geom_line() + 
  facet_wrap(~group, ncol=4) +
  labs(x='White Ethnocentrism', y='Pr(Vote Trump)') +
  scale_x_continuous(limits=c(-0.05,1.05),
                     breaks=c(0,0.25,0.5,0.75,1),
                     labels=c('0','0.25','0.5','0.75','1')) +
  scale_y_continuous(breaks=c(0,0.2,0.4,0.6,0.8,1),
                     limits=c(-0.05,1)) +
  goldenScatterCAtheme +
  geom_point(data=labelling, aes(x=x, y=pe,  group=group), size=7.5, shape=21, fill='white') +
  geom_text(data=labelling, aes(x=x, y=pe,  group=group, label=pe*100),size=2.5) +
  theme(axis.text.y = element_blank())
ggsave(width=8, height=4,filename = 'figures/figure4.eps')

######################
# Figure 5 - Sexism and Trump Support.
######################

set.seed(12345)
form1 <- trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5
form2 <- trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5
form3 <- trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + korean + aapi_other
form4 <- trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other

summary(m1 <- glm(form1, dat1 , family=binomial(link='logit')))
summary(m2 <- glm(form2, dat2, family=binomial(link='logit')))
summary(m3 <- glm(form3, dat3, family=binomial(link='logit')))
summary(m4 <- glm(form4, dat4, family=binomial(link='logit')))

pe1 <- coef(m1); vc1 <- vcov(m1); sb1 <- MASS::mvrnorm(10000, pe1, vc1)
pe2 <- coef(m2); vc2 <- vcov(m2); sb2 <- MASS::mvrnorm(10000, pe2, vc2)
pe3 <- coef(m3); vc3 <- vcov(m3); sb3 <- MASS::mvrnorm(10000, pe3, vc3)
pe4 <- coef(m4); vc4 <- vcov(m4); sb4 <- MASS::mvrnorm(10000, pe4, vc4)

xhyp <- seq(0,1,length.out=100)
mm1 <- cfMake(form1, dat1, nscen=length(xhyp))
mm2 <- cfMake(form2, dat2, nscen=length(xhyp))
mm3 <- cfMake(form3, dat3, nscen=length(xhyp))
mm4 <- cfMake(form4, dat4, nscen=length(xhyp))

mm1 <- cfChange(mm1, 'sexism', x=xhyp, scen=1:length(xhyp))
mm2 <- cfChange(mm2, 'sexism', x=xhyp, scen=1:length(xhyp))
mm3 <- cfChange(mm3, 'sexism', x=xhyp, scen=1:length(xhyp))
mm4 <- cfChange(mm4, 'sexism', x=xhyp, scen=1:length(xhyp))

res1 <- as_tibble(logitsimev(mm1, sb1)) %>% mutate(x=xhyp, group='White')
res2 <- as_tibble(logitsimev(mm2, sb2)) %>% mutate(x=xhyp, group='Black')
res3 <- as_tibble(logitsimev(mm3, sb3)) %>% mutate(x=xhyp, group='Asian')
res4 <- as_tibble(logitsimev(mm4, sb4)) %>% mutate(x=xhyp, group='Latino')

dat <- bind_rows(res1, res2, res3, res4) %>%
  mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))
diff_sex <- dat %>%
  filter(x %in% c(0,1)) %>%
  select(pe, group) %>%
  mutate(z=rep(c('min','max'),4)) %>%
  pivot_wider(id_cols = group,names_from=z,values_from=pe) %>%
  mutate(diff=max-min, outcome='Sexism') %>%
  select(group, diff, outcome)

labelling <- dat %>% filter(x==0 | x==1) %>% mutate(pe = round(pe, 2)) %>%
  mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))
dat %>%
  ggplot(aes(x=x, y=pe, ymin=lower, ymax=upper)) +
  geom_ribbon(alpha=0.3, fill='purple') +
  geom_line() + 
  facet_wrap(~group, ncol=4) +
  labs(x='Sexism', y='Pr(Vote Trump)') +
  scale_x_continuous(limits=c(-0.05,1.05),
                     breaks=c(0,0.25,0.5,0.75,1),
                     labels=c('0','0.25','0.5','0.75','1')) +
  scale_y_continuous(breaks=c(0,0.2,0.4,0.6,0.8,1),
                     limits=c(-0.05,1)) +
  goldenScatterCAtheme +
  geom_point(data=labelling, aes(x=x, y=pe,  group=group), size=7.5, shape=21, fill='white') +
  geom_text(data=labelling, aes(x=x, y=pe,  group=group, label=pe*100),size=2.5) +
  theme(axis.text.y = element_blank())
ggsave(width=8, height=4,filename = 'figures/figure5.eps')

######################
# Appendix Table C2 - Trump Vote:White Ethnocentrism and Sexism
# this code just saves the regressions from above in table format for the appendix
######################

summary(m9 <- glm(trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5, data=df[df$white==1,],family=binomial(link='logit')))
summary(m10 <- glm(trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5, data=df[df$black==1,],family=binomial(link='logit')))
summary(m11 <- glm(trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean, df[df$aapi==1,],family=binomial(link='logit')))
summary(m12 <- glm(trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other, data=df[df$latino==1,],family=binomial(link='logit')))
summary(m13 <- glm(trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5, df[df$white==1,],family=binomial(link='logit')))
summary(m14 <- glm(trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5, df[df$black==1,],family=binomial(link='logit')))
summary(m15 <- glm(trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean, df[df$aapi==1,],family=binomial(link='logit')))
summary(m16 <- glm(trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other, df[df$latino==1,],family=binomial(link='logit')))

stargazer(m9, m10, m11, m12, m13, m14, m15, m16,
          out='tables/c2.tex',
          type="text",
          dep.var.labels = "",
          column.labels = c("White","Black", "AAPI","Latino", 
                            "White","Black", "AAPI","Latino", 
                            "White","Black", "AAPI","Latino"),
          model.names = F,
          omit.stat = c('ll','aic'),
          title = "Logistic Regression Models Predicting Trump Vote",
          covariate.labels = c("White Ethnocentrism","Sexism","College Education",
                               "Household Income","Gender", "Age", "Party ID",
                               "Indian", "Other AAPI", "Korean",
                               "Ideology", "Cuban", "Puerto Rican", "Hispanic Other", 
                              "Constant"))

######################
# Appendix C3 and C4 - Trump Favorability
######################

#Immigration policy + racial resentment
summary(m1 <- lm(trump_fav ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5, data=df[df$white==1,]))
summary(m2 <- lm(trump_fav ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5, data=df[df$black==1,]))
summary(m3 <- lm(trump_fav ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean, df[df$aapi==1,]))
summary(m4 <- lm(trump_fav ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5+ cuban + puertorican + hispanic_other, data=df[df$latino==1,]))
summary(m5 <- lm(trump_fav ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5, data=df[df$white==1,]))
summary(m6 <- lm(trump_fav ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5, data=df[df$black==1,]))
summary(m7 <- lm(trump_fav ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean, df[df$aapi==1,]))
summary(m8 <- lm(trump_fav ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5+cuban + puertorican + hispanic_other, data=df[df$latino==1,]))

stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          type="text",
          out = 'tables/c3.tex',
          dep.var.labels = "",
          column.labels = c("White","Black", "AAPI","Latino", 
                            "White","Black", "AAPI","Latino", 
                            "White","Black", "AAPI","Latino"),
          model.names = F,
          omit.stat = c("ser","f",'rsq'),
          title = "OLS Regression Models Predicting Trump Approval",
          covariate.labels = c("Immigration Policy","Racial Resentment",
                               "College Education","Household Income",
                               "Gender", "Age", "Party ID", "Ideology", 
                               "Indian", "Other AAPI", "Korean",
                               "Cuban", "Puerto Rican", "Hispanic Other", 
                                "Constant"),
          se = starprep(m1, m2, m3, m4, m5, m6, m7, m8))

summary(m9 <- lm(trump_fav ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5, data=df[df$white==1,]))
summary(m10 <- lm(trump_fav ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5, data=df[df$black==1,]))
summary(m11 <- lm(trump_fav ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean, df[df$aapi==1,]))
summary(m12 <- lm(trump_fav ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other, data=df[df$latino==1,]))
summary(m13 <- lm(trump_fav ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5, df[df$white==1,]))
summary(m14 <- lm(trump_fav ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5, df[df$black==1,]))
summary(m15 <- lm(trump_fav ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean, df[df$aapi==1,]))
summary(m16 <- lm(trump_fav ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other, df[df$latino==1,]))

stargazer(m9, m10, m11, m12, m13, m14, m15, m16,
          type="text",
          dep.var.labels = "",
          out='tables/c4.tex',
          column.labels = c("White","Black", "AAPI","Latino",
                            "White","Black", "AAPI","Latino", 
                            "White","Black", "AAPI","Latino"),
          model.names = F,
          omit.stat = c("ser","f",'rsq'),
          title = "OLS Regression Models Predicting Trump Favorability",
          covariate.labels = c("White Ethnocentrism","Sexism","College Education",
                               "Household Income","Gender", "Age", "Party ID","Ideology", 
                               "Indian", "Other AAPI", "Korean",
                               "Cuban", "Puerto Rican", "Hispanic Other", 
                                "Constant"),
          se = starprep(m9, m10, m11, m12, m13, m14, m15, m16))

######################
# Appendix C5 - Trump Vote Model with all IVs
######################

summary(mAll <- glm(trump_hr ~ imm_policy + rr + ethnocentrism2 + sexism + 
                      college_educ + household_income + female + age + 
                      pid7_legacy + ideo5 + white +black + cuban + 
                      puertorican + hispanic_other + indian + aapi_other + 
                      korean, df, family = binomial(link='logit')))

stargazer(mAll,
          type="text",
          out='tables/c5.tex',
          dep.var.labels = "",
          omit.stat=c('aic','ll'),
          model.names = F,
          title = "Logistic Regression Model Predicting Trump Favorability",
          covariate.labels = c("Immigration Policy","Racial Resentment","White Ethnocentrism","Sexism",
                               "College Education","Household Income","Gender", "Age", "Party ID", 
                               "Ideology", "White", "Black", "Cuban", "Puerto Rican", "Hispanic Other", 
                               "Indian", "Other AAPI", "Korean", "Constant"))

######################
# Figure 6 - Attention, Status Hierarchy, and Support for Trump.
######################

goldenScatterCAthemeFlip <- theme(
  ## Removes main plot gray background
  panel.background = element_rect(fill = "white"),
  ## Golden rectangle plotting area (leave out for square)
  aspect.ratio = ((1 + sqrt(5))/2)^(-1),
  ## All axes changes
  axis.ticks.length = unit(0, "char"), # longer ticks
  ## Horizontal axis changes
  axis.text.y = element_text(color = "black", size = 12),
  ## match type of axis labels and titles
  axis.title.y = element_text(size = 12,
                              margin = margin(t = 7.5, r = 0, b = 0, l = 0)),
  ## match type; pad space between title and labels
  ## Vertical axis changes
  axis.ticks.y = element_blank(), # no y axis ticks (gridlines suffice)
  axis.text.x = element_text(color = "black", size = 12,
                             margin = margin(t = 0, r = -4, b = 0, l = 0)),
  ## match type of axis labels and titles, pad
  axis.title.x = element_text(size = 12,
                              margin = margin(t = 0, r = 7.5, b = 0, l = 0)),
  ## match type of axis labels and titles, pad
  ## Legend
  legend.key = element_rect(fill = NA, color = NA),
  ## Remove unhelpful gray background
  ## Gridlines (in this case, horizontal from left axis only
  panel.grid.major.x = element_line(color = "gray45", size = 0.2),
  panel.grid.major.y = element_blank(),
  ## Faceting (small multiples)
  strip.background = element_blank(),
  ## Remove unhelpful trellis-like shading of titles
  strip.text.x = element_text(size=12), # Larger facet titles
  strip.text.y = element_blank(), # No titles for rows of facets
  strip.placement = "outside", # Place titles outside plot
  panel.spacing.x = unit(1.25, "lines"), # Horizontal space b/w plots
  panel.spacing.y = unit(1, "lines") # Vertical space b/w plots
)
class(goldenScatterCAthemeFlip) 

# interest interaction

dat1 <- df[df$white==1,]
dat2 <- df[df$black==1,]
dat3 <-  df[df$aapi==1,]
dat4 <- df[df$latino==1,]

diffPlotInterest <- function(form1, form2, form3, form4,outcome ,var){

  summary(m1 <- glm(form1, dat1, family=binomial))
  summary(m2 <- glm(form2, dat2, family=binomial))
  summary(m3 <- glm(form3, dat3, family=binomial))
  summary(m4 <- glm(form4, dat4, family=binomial))
  
  pe1 <- coef(m1); vc1 <- vcov(m1); sb1 <- MASS::mvrnorm(10000, pe1, vc1)
  pe2 <- coef(m2); vc2 <- vcov(m2); sb2 <- MASS::mvrnorm(10000, pe2, vc2)
  pe3 <- coef(m3); vc3 <- vcov(m3); sb3 <- MASS::mvrnorm(10000, pe3, vc3)
  pe4 <- coef(m4); vc4 <- vcov(m4); sb4 <- MASS::mvrnorm(10000, pe4, vc4)
  
  mm1 <- cfMake(form1, dat1, nscen=2)
  mm2 <- cfMake(form2, dat2, nscen=2)
  mm3 <- cfMake(form3, dat3, nscen=2)
  mm4 <- cfMake(form4, dat4, nscen=2)
  
  mm1 <- cfChange(mm1, var, xpre=0, x=1, scen=1:2)
  mm2 <- cfChange(mm2, var, xpre=0, x=1, scen=1:2)
  mm3 <- cfChange(mm3, var, xpre=0, x=1, scen=1:2)
  mm4 <- cfChange(mm4, var, xpre=0, x=1, scen=1:2)
  
  mm1 <- cfChange(mm1, 'attention', xpre=0, x=0, scen=1)
  mm2 <- cfChange(mm2, 'attention', xpre=0, x=0, scen=1)
  mm3 <- cfChange(mm3, 'attention', xpre=0, x=0, scen=1)
  mm4 <- cfChange(mm4, 'attention', xpre=0, x=0, scen=1)
  
  mm1 <- cfChange(mm1, 'attention', xpre=1, x=1, scen=2)
  mm2 <- cfChange(mm2, 'attention', xpre=1, x=1, scen=2)
  mm3 <- cfChange(mm3, 'attention', xpre=1, x=1, scen=2)
  mm4 <- cfChange(mm4, 'attention', xpre=1, x=1, scen=2)
  
  res1 <- as_tibble(logitsimfd(mm1, sb1)) %>% 
    mutate(group='White', outcome=outcome, Attention = c("Low",'High'))
  res2 <- as_tibble(logitsimfd(mm2, sb2)) %>% 
    mutate(group='Black',outcome=outcome, Attention = c("Low",'High'))
  res3 <- as_tibble(logitsimfd(mm3, sb3)) %>% 
    mutate(group='Asian', outcome=outcome, Attention = c("Low",'High'))
  res4 <- as_tibble(logitsimfd(mm4, sb4)) %>% 
    mutate(group='Latino', outcome=outcome, Attention = c("Low",'High'))
  
  dat <- bind_rows(res1, res2, res3, res4) %>%
    dplyr::rename(diff=pe) %>%
    mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))
  return(dat)
}

form1 <- trump_hr ~ imm_policy*attention + college_educ + household_income + female + age + pid7_legacy + ideo5
form2 <- trump_hr ~ imm_policy*attention + college_educ + household_income + female + age + pid7_legacy + ideo5
form3 <- trump_hr ~ imm_policy*attention + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form4 <- trump_hr ~ imm_policy*attention + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form5 <- trump_hr ~ sexism*attention + college_educ + household_income + female + age + pid7_legacy + ideo5
form6 <- trump_hr ~ sexism*attention + college_educ + household_income + female + age + pid7_legacy + ideo5
form7 <- trump_hr ~ sexism*attention + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form8 <- trump_hr ~ sexism*attention + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form9 <- trump_hr ~ rr*attention + college_educ + household_income + female + age + pid7_legacy + ideo5
form10 <- trump_hr ~ rr*attention + college_educ + household_income + female + age + pid7_legacy + ideo5
form11 <- trump_hr ~ rr*attention + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form12 <- trump_hr ~ rr*attention + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form13 <- trump_hr ~ ethnocentrism2*attention + college_educ + household_income + female + age + pid7_legacy + ideo5
form14 <- trump_hr ~ ethnocentrism2*attention + college_educ + household_income + female + age + pid7_legacy + ideo5
form15 <- trump_hr ~ ethnocentrism2*attention + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form16 <- trump_hr ~ ethnocentrism2*attention + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other

summary(m1 <- glm(form1, dat1, family=binomial))
summary(m2 <- glm(form2, dat2, family=binomial))
summary(m3 <- glm(form3, dat3, family=binomial))
summary(m4 <- glm(form4, dat4, family=binomial))
summary(m5 <- glm(form5, dat1, family=binomial))
summary(m6 <- glm(form6, dat2, family=binomial))
summary(m7 <- glm(form7, dat3, family=binomial))
summary(m8 <- glm(form8, dat4, family=binomial))
summary(m9 <- glm(form9, dat1, family=binomial))
summary(m10 <- glm(form10, dat2, family=binomial))
summary(m11 <- glm(form11, dat3, family=binomial))
summary(m12 <- glm(form12, dat4, family=binomial))
summary(m13 <- glm(form13, dat1, family=binomial))
summary(m14 <- glm(form14, dat2, family=binomial))
summary(m15 <- glm(form15, dat3, family=binomial))
summary(m16 <- glm(form16, dat4, family=binomial))

interest <- bind_rows(
  diffPlotInterest(form1, form2, form3, form4,'Immigration', 'imm_policy'),
  diffPlotInterest(form5, form6, form7, form8,'Sexism', 'sexism'),
  diffPlotInterest(form9, form10, form11, form12,'Racial Resentment', 'rr'),
  diffPlotInterest(form13, form14, form15, form16,'White Ethnocentrism', 'ethnocentrism2')
) 

ggplot(interest, aes(x=reorder(outcome,diff), y=diff, ymin=lower, ymax=upper, shape=Attention)) +
  geom_errorbar(width=0, position=position_dodge(width=0.2)) +
  geom_point(size=2, fill='white', position=position_dodge(width=0.2)) + 
  goldenScatterCAthemeFlip +
  facet_wrap(~group,ncol=2) +
  coord_flip() +
  labs(x='', y='Change Pr(Vote Trump)')  +
  scale_shape_manual(values=c(21,24), name='Attention') +
  theme(legend.position='bottom') +
  scale_y_continuous(labels=scales::percent, limits=c(0,1))
ggsave(width=6, height=5,filename = 'figures/figure6.eps')

# Appendix Tables C11 and C12
stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          out = 'tables/c11.tex',
          type="text",
          omit.stat = c('ll','aic'),
          dep.var.labels = "",
          column.labels = c("White","Black", "AAPI", "Latino","White","Black", "AAPI", "Latino"),
          model.names = F,
          title = "Attention, Status Hierarchy, and Support for Trump: Immigration & Sexism",
          covariate.labels = c("Immigration Policy","Sexism", "Attention", "College Education",
                               "Household Income","Gender", "Age", "Party ID", "Ideology", 
                               "Indian", "Other AAPI", "Korean", 
                               "Cuban", "Puerto Rican", "Hispanic Other", 
                               "Immigration Policy*Attention","Sexism*Attention", "Constant"))
stargazer(m9, m10, m11, m12, m13, m14, m15, m16,
          type="text",
          out = 'tables/c12.tex',
          omit.stat = c('ll','aic'),
          dep.var.labels = "",
          column.labels = c("White","Black", "AAPI", "Latino",
                            "White","Black", "AAPI", "Latino"),
          model.names = F,
          title = "Attention, Status Hierarchy, and Support for Trump: Racial Resentment & White Ethnocentrism",
          covariate.labels = c("Racial Resentment","White Ethnocentrism", "Attention","College Education",
                               "Household Income","Gender", "Age", "Party ID", "Ideology", 
                               "Indian", "Other AAPI", "Korean",
                               "Cuban", "Puerto Rican", "Hispanic Other",
                               "Racial Resentment*Attention","White Ethnocentrism*Attention", "Constant"))

######################
# Figure 7 and Appendix Tables C13 - C16 - Status Hierarchy and Support for Other Republicans.
######################

diffPlot <- function(form1, form2, form3, form4,outcome ,var){
  set.seed(12345)
  summary(m1 <- glm(form1, dat1, family=binomial))
  summary(m2 <- glm(form2, dat2, family=binomial ))
  summary(m3 <- glm(form3, dat3, family=binomial ))
  summary(m4 <- glm(form4, dat4, family=binomial ))
  pe1 <- coef(m1); vc1 <- vcov(m1)
  sb1 <- MASS::mvrnorm(10000, pe1, vc1)
  pe2 <- coef(m2); vc2 <- vcov(m2)
  sb2 <- MASS::mvrnorm(10000, pe2, vc2)
  pe3 <- coef(m3); vc3 <- vcov(m3)
  sb3 <- MASS::mvrnorm(10000, pe3, vc3)
  pe4 <- coef(m4); vc4 <- vcov(m4)
  sb4 <- MASS::mvrnorm(10000, pe4, vc4)
  mm1 <- cfMake(form1, dat1, nscen=1)
  mm2 <- cfMake(form2, dat2, nscen=1)
  mm3 <- cfMake(form3, dat3, nscen=1)
  mm4 <- cfMake(form4, dat4, nscen=1)
  mm1 <- cfChange(mm1, var, xpre=0, x=1, scen=1)
  mm2 <- cfChange(mm2, var, xpre=0, x=1, scen=1)
  mm3 <- cfChange(mm3, var, xpre=0, x=1, scen=1)
  mm4 <- cfChange(mm4, var, xpre=0, x=1, scen=1)
  res1 <- as_tibble(logitsimfd(mm1, sb1)) %>% mutate(group='White')
  res2 <- as_tibble(logitsimfd(mm2, sb2)) %>% mutate(group='Black')
  res3 <- as_tibble(logitsimfd(mm3, sb3)) %>% mutate(group='Asian')
  res4 <- as_tibble(logitsimfd(mm4, sb4)) %>% mutate(group='Latino')
  dat <- bind_rows(res1, res2, res3, res4) %>%
    dplyr::rename(diff=pe) %>%
    mutate(outcome=outcome) %>%
    mutate(group=factor(group, levels=c('White','Black','Latino','Asian')))
  return(dat)
}

form1 <- trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5
form2 <- trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5
form3 <- trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + chinese + korean
form4 <- trump_hr ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form5 <- trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5
form6 <- trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5
form7 <- trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + chinese + korean
form8 <- trump_hr ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form9 <- trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5
form10 <- trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5
form11 <- trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + chinese + korean
form12 <- trump_hr ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form13 <- trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5
form14 <- trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5
form15 <- trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + chinese + korean
form16 <- trump_hr ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other

trump_vote_pe <- bind_rows(
  diffPlot(form1, form2, form3, form4,'Immigration', 'imm_policy'),
  diffPlot(form5, form6, form7, form8,'Sexism', 'sexism'),
  diffPlot(form9, form10, form11, form12,'Racial Resentment', 'rr'),
  diffPlot(form13, form14, form15, form16,'White Ethnocentrism', 'ethnocentrism2')
) %>%
  mutate(z='Trump Vote')


form1 <- vote_romney ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5
form2 <- vote_romney ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5
form3 <- vote_romney ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form4 <- vote_romney ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form5 <- vote_romney ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5
form6 <- vote_romney ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5
form7 <- vote_romney ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form8 <- vote_romney ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form9 <- vote_romney ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5
form10 <- vote_romney ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5
form11 <- vote_romney ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form12 <- vote_romney ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form13 <- vote_romney ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5
form14 <- vote_romney ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5
form15 <- vote_romney ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form16 <- vote_romney ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other

# C13
summary(m1 <- glm(form1, dat1, family=binomial))
summary(m2 <- glm(form2, dat2, family=binomial))
summary(m3 <- glm(form3, dat3, family=binomial))
summary(m4 <- glm(form4, dat4, family=binomial))
summary(m5 <- glm(form5, dat1, family=binomial))
summary(m6 <- glm(form6, dat2, family=binomial))
summary(m7 <- glm(form7, dat3, family=binomial))
summary(m8 <- glm(form8, dat4, family=binomial))

stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          type="text",
          dep.var.labels = "",
          omit.stat = c('ll','aic'),
          out='tables/c13.tex',
          column.labels = c("White","Black", "AAPI", "Latino",
                            "White","Black", "AAPI", "Latino"),
          model.names = F,
          title = "Logistic Regression Models Predicting Romney Vote",
          covariate.labels = c("Immigration Policy","Sexism","College Education",
                               "Household Income","Female",
                               "Age", "Party ID", "Ideology",  
                               "Indian", "Other AAPI", "Korean",
                               "Cuban", "Puerto Rican", "Hispanic Other",
                               "Constant"))
# C14
summary(m9 <-  glm(form9, dat1, family=binomial))
summary(m10 <- glm(form10, dat2, family=binomial))
summary(m11 <- glm(form11, dat3, family=binomial))
summary(m12 <- glm(form12, dat4, family=binomial))
summary(m13 <- glm(form13, dat1, family=binomial))
summary(m14 <- glm(form14, dat2, family=binomial))
summary(m15 <- glm(form15, dat3, family=binomial))
summary(m16 <- glm(form16, dat4, family=binomial))

stargazer(m9, m10, m11, m12, m13, m14, m15, m16,
          type="text",
          dep.var.labels = "",
          omit.stat = c('ll','aic'),
          out='tables/c14.tex',
          column.labels = c("White","Black", "AAPI", "Latino",
                            "White","Black", "AAPI", "Latino"),
          model.names = F,
          title = "Logistic Regression Models Predicting Romney Vote",
          covariate.labels = c("Racial Resentment","White Ethnocentrism","College Education","Household Income","Female",
                               "Age", "Party ID", "Ideology",  
                               "Indian", "Other AAPI", "Korean",
                               "Cuban", "Puerto Rican", "Hispanic Other",
                               "Constant"))

romney_vote_pe <- bind_rows(
  diffPlot(form1, form2, form3, form4,'Immigration', 'imm_policy'),
  diffPlot(form5, form6, form7, form8,'Sexism', 'sexism'),
  diffPlot(form9, form10, form11, form12,'Racial Resentment', 'rr'),
  diffPlot(form13, form14, form15, form16,'White Ethnocentrism', 'ethnocentrism2')
) %>%
  mutate(z='Romney Vote')

form1 <- vote_rep_house ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5
form2 <- vote_rep_house ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5
form3 <- vote_rep_house ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form4 <- vote_rep_house ~ imm_policy + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form5 <- vote_rep_house ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5
form6 <- vote_rep_house ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5
form7 <- vote_rep_house ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form8 <- vote_rep_house ~ sexism + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form9 <- vote_rep_house ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5
form10 <- vote_rep_house ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5
form11 <- vote_rep_house ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form12 <- vote_rep_house ~ rr + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other
form13 <- vote_rep_house ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5
form14 <- vote_rep_house ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5
form15 <- vote_rep_house ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + indian + aapi_other + korean
form16 <- vote_rep_house ~ ethnocentrism2 + college_educ + household_income + female + age + pid7_legacy + ideo5 + cuban + puertorican + hispanic_other

summary(m1 <- glm(form1, dat1, family=binomial))
summary(m2 <- glm(form2, dat2, family=binomial))
summary(m3 <- glm(form3, dat3, family=binomial))
summary(m4 <- glm(form4, dat4, family=binomial))
summary(m5 <- glm(form5, dat1, family=binomial))
summary(m6 <- glm(form6, dat2, family=binomial))
summary(m7 <- glm(form7, dat3, family=binomial))
summary(m8 <- glm(form8, dat4, family=binomial))

# C15
stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          type="text",
          dep.var.labels = "",
          omit.stat = c('ll','aic'),
          out='tables/c15.tex',
          column.labels = c("White","Black", "AAPI", "Latino",
                            "White","Black", "AAPI", "Latino"),
          model.names = F,
          title = "Logistic Regression Models Predicting House Republican Vote",
          covariate.labels = c("Immigration Policy","Sexism","College Education",
                               "Household Income","Gender", "Age", "Party ID", "Ideology", 
                               "Indian", "Other AAPI", "Korean",
                               "Cuban", "Puerto Rican", "Hispanic Other",
                               "Constant"))

# C16
summary(m9 <-  glm(form9, dat1, family=binomial))
summary(m10 <- glm(form10, dat2, family=binomial))
summary(m11 <- glm(form11, dat3, family=binomial))
summary(m12 <- glm(form12, dat4, family=binomial))
summary(m13 <- glm(form13, dat1, family=binomial))
summary(m14 <- glm(form14, dat2, family=binomial))
summary(m15 <- glm(form15, dat3, family=binomial))
summary(m16 <- glm(form16, dat4, family=binomial))

stargazer(m9, m10, m11, m12, m13, m14, m15, m16,
          type="text",
          dep.var.labels = "",
          omit.stat = c('ll','aic'),
          out='tables/c16.tex',
          column.labels = c("White","Black", "AAPI", "Latino",
                            "White","Black", "AAPI", "Latino"),
          model.names = F,
          title = "Logistic Regression Models Predicting House Republican Vote",
          covariate.labels = c("Racial Resentment","White Ethnocentrism","College Education",
                               "Household Income","Gender", "Age", "Party ID", "Ideology", 
                               "Indian", "Other AAPI", "Korean",
                               "Cuban", "Puerto Rican", "Hispanic Other",
                               "Constant"))

house_vote_pe <- bind_rows(
  diffPlot(form1, form2, form3, form4,'Immigration', 'imm_policy'),
  diffPlot(form1=form5, form2=form6, form3=form7, form4=form8,'Sexism', 'sexism'),
  diffPlot(form9, form10, form11, form12,'Racial Resentment', 'rr'),
  diffPlot(form13, form14, form15, form16,'White Ethnocentrism', 'ethnocentrism2')
) %>%
  mutate(z='House Vote Rep')

bind_rows(trump_vote_pe, romney_vote_pe, house_vote_pe) %>%
  mutate(z = factor(z, levels=c(
    'Romney Vote','House Vote Rep','Trump Vote'
  ))) %>%
  ggplot(aes(x=reorder(outcome,diff), y=diff, ymin=lower, ymax=upper, shape=z)) +
  geom_errorbar(width=0, position=position_dodge(width=0.2)) +
  geom_point(size=2, fill='white', position=position_dodge(width=0.2)) + 
  goldenScatterCAthemeFlip +
  facet_wrap(~group,ncol=2) +
  coord_flip() +
  labs(x='', y='Change Pr(Vote Trump)')  +
  scale_shape_manual(values=c(21,24,19), name='') +
  theme(legend.position='bottom') +
  scale_y_continuous(labels=scales::percent)
ggsave(width=6, height=5,filename = 'figures/figure7.eps')

